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We study the existence, stability, and dynamics of symmetric and anti-symmetric states of quasi- 
one-dimensional polariton condensates in double-well potentials, in the presence of nonresonant 
pumping and nonlinear damping. Some prototypical features of the system, such as the bifurcation 
of asymmetric solutions, are similar to the Hamiltonian analog of the double- well system considered 
in the realm of atomic condensates. Nevertheless, there are also some nontrivial differences including, 
e.g., the unstable nature of both the parent and the daughter branch emerging in the relevant 
pitchfork bifurcation for slightly larger values of atom numbers. Another interesting feature that 
does not appear in the atomic condensate case is that the bifurcation for attractive interactions 
is slightly sub-critical instead of supercritical. These conclusions of the bifurcation analysis are 
corroborated by direct numerical simulations examining the dynamics of the system in the unstable 
regime. 

I. INTRODUCTION 

Over the past few years, a novel direction in the study of Bose-Einstein condensation has captured a considerable 
amount of attention. This concerns the observation of exciton-polariton Bose-Einstein condensates (BECs) in semi- 
conductor microcavities [H-Q- A fundamental feature of these exciton-polariton BECs is that, upon confinement, 
the excitons (bound pairs of electrons and holes) couple strongly to the incident light creating the polariton quasi- 
particles @ . The resulting exciton-polariton BEC possesses a number of remarkable properties that we briefly touch 
upon below. 

The radiative lifetime of the polaritons is the shorter relaxation time scale of the system being of the order of 
1-10 ps @. On the other hand, the light mass of the exciton-polaritons provides this system with a significantly 
higher condensation temperature. The photonic component of the exciton-polaritons is responsible for their short 
lifetime which, in turn, does not allow thermalization; instead, it produces a non-equilibrium condensate, wherein the 
presence of external pumping from an exciton reservoir is critical towards a counter-balance of the polariton loss. In 
such genuinely non-equilibrium condensates, numerous remarkable features have been not only theoretically predicted 
but also experimentally established; these include the fiow without scattering (analog of the flow without friction) Q, 
the existence of vortices [1] (see also Ref. @ for vortex dipole dynamics and Ref. fiol for observations thereof), the 
collective dynamics [ll|, as well as remarkable applications such as spin switches |l2| and light emitting diodes [l^ 
operating even near room temperatures. 

Perhaps the most customary approach to modeling exciton-polariton BECs involves the coupling of the evolution 
of the polaritons to that of the exciton reservoir which enables their production (and which features diffusive spatial 
dynamics of the excitons); this way, the model takes the form of two coupled complex Ginzburg-Landau (cGL) 
equations describing the evolution of exciton and photon wavcfunctions [lj-|l6[. Nevertheless, it has been proposed 
m Refs. [IMl that a single cGL equation for the macroscopically occupied polariton state can also be used in 
a way consistent with experimental observations [1^. A similar approach was followed in Ref. [2l| where a BEC 
of magnon quasi-particles, incorporating a source term rather than an amplification of the field, was shown to be 
phenomenologically described by a system two nonlinearly coupled cGL-type equations. In the context of the single 
cGL model for the polaritons, there exists a localized (pumping) region of gain and a nonlinear saturating loss 
term, in addition to all the standard terms (quantum pressure, external parabolic trapping and repulsive interatomic 
interaction) that one encounters in atomic BECs [13 . Furthermore, it should be pointed out that the prototypical 
setting where experiments have been conducted is two-dimensional in nature. Yet, highly anisotropic traps (similar 
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to what has been done in atomic BECs [2^) can be envisioned which reduce the effective dynamics to a quasi one- 
dimensional (ID) setting [23l - [28| . Moreover, recent experimental advances have enabled the use of thin microwires 
in order to guide the condensates along the direction of the wire [l^- In this setting, the recent analysis of Ref. [s^l 
presented a number of striking characteristics due to the interplay of gain and loss terms with the standard ones of 
atomic BECs. Prominent examples included the destabilization of the nodeless state of the system and the creation 
of stability inversions (where states with nodes would be more robust), as well as the existence of bubble-like and 
sawtooth-like solutions in the system. 

A very interesting research direction in the physics of atomic and polariton BECs concerns the dynamics of the 
condensates in a double- well potential. The latter can be created in atomic BEC experiments through the combination 
of a parabolic trap and a p eriodic (so-called optical lattice) potenti al g enerated through the interference of laser beams 
illuminating the BEC [3l[. Relevant experiments in atomic BECs j32l [33| have paved the way towards the exploration 
of numerous features such as tunneling and Josephson oscillations for small numbers of atoms in the condensate, and 
macroscopic quantum self-trapped states, as well as symmetry-breaking effects for large atom numbers. On the other 
hand, double-well potentials can also be created in polariton BEC experiments in microcavities by applying stress 
1l 2JI J by employing photolithographic techniques (23l . [2^ , or allowing natural formation during the sample growth 
35 1 . Importantly, the latter technique was used for the study of a "polariton Josephson junction" [35j . in the spirit of 
earlier studies on "bosonic Josephson junctions" [s^ in the context of atomic BECs. Importantly, a large volume of 
theoretical studies has accompanied these developments, first in the context of atomic BECs, through investigations 
related to finite-mode reductions and symmetry-breaking bifurcations [37i - |4^ . quantum effects [45ij . and nonlinear 
variants of the double-well potential 1461 . and more recently in the context of polariton condensates, especially as 
concerns Josephson oscillations therein [47[. It should be mentioned in passing that similar (spontaneous symmetry 
breaking) effects have been monitored in the realm of nonlinear optics: in this context, formation of asymmetric 
states in dual-core fibers [48| , self-guided laser beams in Kerr media j49| , and optically-induced dual-core waveguiding 
structures in photorefractive crystals [soj have been reported. 

It is the aim of the present work to combine these two themes, namely the focus on the cxciton-polariton BEC 
with pumping and loss and the fundamental interest in the understanding of double-well trapping potentials in a 
spirit similar to the proposal of Ref. [iJI. In particular, we will consider the single-component model of Refs. [l7l - [l9| 
combined with a double- well potential in a quasi-lD (e.g., microwire) setting. We will attempt a systematic (Galerkin) 
finite-dimensional reduction of the system via projection to the two principal cigenstatcs of the potential, and will 
derive a damped-driven system of ordinary differential equations (ODEs) that have been shown in the Hamiltonian 



case to capture the essence of the statics [51[ and dynamics [52| of double-well potentials. We will then examine 
the bifurcation structure of the resulting ODEs and compare it to that of the original partial differential equation 
(PDE) model. This already provides us with a number of interesting features that distinguish this system from 
its Hamiltonian analog. For instance, in the case of attractive interatomic interactions (which is studied together 
with that of repulsive interactions) the relevant symmetry-breaking pitchfork bifurcation is subcritical instead of 
supercritical as in the Hamiltonian case. Furthermore, both branches that emerge from the pitchfork bifurcations, the 
stable asymmetric one and the (now) unstable "parent" branch, both appear to become destabilized in this polariton 
BEC setting for slightly larger nonlinearitics, posing the natural question of what is the stable dynamics for larger 
values of the nonlinearity. These questions will in part be addressed via direct numerical simulations. 

Our presentation will be structured as follows. First, in Section II, we will present the model and its theoretical 
study via the Galerkin analysis. In Section HI, we will study the model numerically and compare the results of the 
numerical bifurcation analysis with the prediction of the Galerkin approximation. We will also complement these 
results with direct numerical simulations of the original model. Finally, in section IV, we summarize our results and 
present our conclusions. 



II. MODEL SETUP AND ANALYTICAL PREDICTIONS 

In our analysis below, we adopt the model of Refs. [T7l - [l9j . It has been argued in these works that the original 
exciton-polariton system given by a set of two coupled equations can be effectively reduced to a single cGL equation 
with a nonlinear saturating loss term. This reduction can be used when the reservoir mean-field potential is negligible 
and the spot size is large compared with the condensate size (i.e., if we can consider that the spot width is the same 
of the spatial extent of the system). In particular, the amplification of the existing field introduces a gain and hence 
acts as a generator of polaritons. Then the loss term saturates this gain beyond a certain threshold. These two terms 
are analogous to the pumping of polaritons from the excitons and to the natural decay of the polaritons. This reduced 
model can be expressed in dimensionless form as follows: 



idtu ~ —d^u + sjupit -I- V{x)u ~ iju + i [x{x) — cr\u\'^'\ u. 



(1) 
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The above model is actually a complex Ginzburg-Landau equation [53| for the complex order parameter u{x, t), which 
is assumed to evolve in the presence of the effectively-lD double- well potential V{x). Equation ([1]) can be applied 
to both the contexts of atomic and polariton BECs: in the first case, the two last terms in the right-hand side of 
Eq. ^ are absent, and the model — known as the Gross-Pitaevskii equation — describes the evolution of the 
macroscopic wavefunction for the cold atoms and ^ is the chemical potential; in the second case, u{x, t) denotes the 
polariton wavefunction, and the last two terms in the right-hand side are included in the model. More specifically, in 
the context of polariton condensates, Eq. ([I]) incorporates (a) the spatially dependent gain term of the form 

X{x) = ae{x,n - \x\), (2) 

where is the step function generating a symmetric pumping spot of "radius" x^ and strength a for the gain, and 
(b) a nonlinear saturation loss term, characterized by its strength a. As concerns the parameter s = ±1, it sets 
the type of nonlinearity (i.e., the type of interactions between atoms or polaritons): for s = +1 the nonlinearity is 
defocusing (i.e., the interactions are repulsive), while for s = — 1 the nonlinearity is focusing (i.e., the interactions are 
attractive). In the context of atomic BECs, the value of s depends on the atom species (e.g., s = -t-1 for ^^Rb or ^■^Na, 
while s = — 1 for ^Li or ®^Rb atoms). On the other hand, in the context of polariton condensates, the sign of the 
effective mass of polaritons [i.e., the sign of the first term in the right-hand side of Eq. ^] may become either positive 
or negative, depending on the values of transverse momentum: in fact, the transition from positive to negative mass 
is associated with the inflection point of the energy-momentum diagram (53 |. Here, we will consider both cases of 
s ~ ±1 to take into regard that the effective polariton mass may be positive or negative, respectively. We finally 
note that the relevant physical time and space scales, as well as physically relevant parameter values associated with 
Eq. H]), can be found in Ref. [T7| . 

In what follows, we will use the Galerkin (few mode truncation) approach of Ref. We start by considering the 
corresponding linear eigenproblem which reads: 

Hu = -dlu + V[x)u = uju, (3) 

whose spectrum consists of a ground state, uq{x), and excited states, Ui{x) (with i > 1). Then, in the weakly nonlinear 
regime, we consider a superposition of the two lowest linear eigenmodes, 

u{x,t) = CQ{t)uo{x) + ci{t)ui{x), (4) 

where co,i(t) are unknown time-dependent complex prefactors; obviously, the above ansatz is relevant for values of 
the chemical potential fi such that higher order modes can be safely ignored. Substituting this ansatz into Eq. ([T]) we 
obtain: 

i {cqUo + dim) = (i^o - m)"oCo + (wi - m)"iCi + s\uf (cqUq + ciui) + i [x{x) - a\u\'^] {cqUq + ciui) , (5) 

where the \u\'^ has not been expanded only for reasons of compactness but should actually be thought as expanded 
according to Eq. (|4]). Next, projecting on uq and ui (i.e., multiplying the above equation by uq and ui and integrating 
over x), and using the orthogonality of the states Ui, we respectively derive the following equations: 

ico = {uJo - fj, + iao)co + [s - ia) {Ao|co|^co -|- (cgC* -I- 2|co|^ci) Tq + (2|ci|^co -I- c^Cq) B + Icil^ciTi} , (6) 

and 

ici = (wi - n + iai)ci + {s - ia) {To\co\^co + {clcl + 2\cof a) B + {2\cifco + clco)Ti + \ci\^ciAi}. (7) 

In the above equations, overdots denote time derivatives, the involved constants (depending on the eigenbasis {ui}) 
take the values Aq = J u^dx, Ai = J ufdx, B = J u^u\dx, Fq = / uiu^dx, and Fi = / uoufdx, while the effective 
gain coefficients read: ao = J x{^)uQdx and ai — J x{x)uidx. We now use amplitude and phase variables for the 
time-dependent prefactors, i.e., Ci = piC^"^' (with the amplitudes pi and phases (f>i being real functions), to derive a 
set of four equations for the unknown functions po,i and 00,1. Introducing the relative phase of the first two modes 
as (/9 = 01 — 00, the above mentioned set of equations takes the following form: 

Pq =aopo - o- {AqpI + 2Bp\pa) + s (FipJ + Taplpi) sin (p 

+ sBp^po sin 2lp — a [Tip'i + 3FoPoPij cos ip — aBp^po cos 2ip, 



00 = - {ujQ ~ p) - s {AqpI + 2Bpf) - a (FopoPi + TiPi/po) sin((5 
— cfBp\ sin 2lp — s (3FojOoPi + TiPi/po) cos — sBp^ cos 2(^3, 



(9) 



Pi =aipi - a {AipI + 2Bplpi) - s {T^pl + Tip\po) sirup 

— sBpqPi siii2(/3 — cr (FoPo + '^^iPiPoj cosip — aBpQpi cos2(^, 

and 

4>i = ~ {uji - p) - s {Aipf + 2Bpl) + cr (ripipo +roPo/pi) sin ((9 
+ cri?/9Q sin2(p — s (SFipipo + roPo/pi) coscp — s-Bpo cos2(^. 

Subtracting Eq. ([9]) from Eq. we can readily obtain an equation for (p, namely: 

ip^^^u-s {A^pI - A^pI) -sB[2 + cos 2^] {pi - p^) - (ropg(p§ - 3p?) + rip2(3p2 „ ^2)-, 

+ ^5^ {TopI{pI + pi) + Tipl{pl+ pI)) + aBsin2^{pl + pi), 
PoPi 



(10) 



(11) 



(12) 



where Aw = wi — a;2. This way, we have arrived to a system of three equations [cf. Eqs. (jH)), ((T0| and ([T2|)] for the 
unknown functions po,i and These equations are subject to an additional constraint stemming from the balance 
condition dN/dt = 0, where N = \u\'^dx is the number of polaritons (mathematically the squared norm). The 
evolution of the latter, can readily be found by multiplying Eq. ([T]) by u*, the complex conjugate of Eq. ([T]) by u, and 
then adding and integrating the resulting equations. It is straightforward to find that the condition for equilibrium 
is: 

{x{x) - <7\u\^)\u\^dx = 0. (13) 

-CJO 

Substituting Eq. ^ into Eq. (|13p , also using the polar decomposition for a (t) [and assuming a definite — even in our 
considerations — parity for the function x(2^)]i find that the balance condition (|13p takes the form: 

(ofoPo + "iPi) - (^oPo + Pi^i + ^pIpIb) - 4cr (pgpiEo + plpo^i) cos ip - 2aplpiB cos 2ip = 0, (14) 

which essentially fixes pi once po and Lp are found and thus reducing the effective number of degrees of freedom for 
our approximations to only two (po and tp). 

Below, we will consider the case of a symmetric double- well potential, for which Fi = Fq = 0. In this case, Eqs. ([U, 
(ITUl) and are reduced to the following simpler form, 

po = aopo - a {AopI + 2Bpipn) + sBpipo sin 2p - crBpipo cos 2ip, (15) 

pi = aipi - a [AipI + 2Bplpi) - sBp^pi sm2p - aBp^pi cos2(p, (16) 

ip = -Auj - s {Aipi - Aapl) - sB[2 + cos2i^] {pi - pi) + crB sin 2<^(p^ + pf), (17) 

while the equilibrium condition is accordingly simplified as: 

(aoPo + aip\) - a {AqpI + p{Ai + ^pIp\B) - 2aplplB cos 2^ = 0. (18) 

We can now turn to the study of stationary solutions (i.e., po ^ pi = ^ 0) resulting from the Galerkin truncation 
analysis. Particularly, from Eq. ((TS|) we obtain two possible solutions: 

) ao-cr {AqpI + 2Bpi) + sBpj sin 2ip - aBpj cos 2ip = 0, 



while from Eq. (1161) we obtain: 



i) Pi = Q 
a 



) ai~ a {Aip\ + 2Bpl) ~ sBplsu\2p- i7Bplcos2ip =0. ^'^^^ 
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Next, multiplying the nontrivial equilibria of Eq. (|T9)) by pg. the one from Eq. ((20|) by pf, and adding the resulting 
equations, we obtain: 



{aopl + aipl) - a {App^ + A^pf + ABplpl) 
2aBplpl 



cos2v3= ^ ^ ^, (21) 



while subtracting Eq. ((20)) from Eq. (jl9p yields: 

a (Aipt - ^oPo + 2B(p2 - ^ ^(^2 _ ^Qg2^) + sB{pI + p^) sin2^ + (ao - ai) = 0. (22) 
Combining now Eq. with Eq. (|17p we finally obtain the result: 



2 2\ ■ r, (jAlu — s{aQ — ai) 
__ 



{pi+Pi)su.2^= "^^^,2:\.r - (23) 



Let us now focus again on Eqs. (|T5|) and (|T6l) : it is clear that if Eq. (fT6)) is satisfied for pi = then = ^3^' ^^^"^ 
Eq. ([TS]) is satisfied with po = then pf = ^-j-- Aside from these trivial symmetric and anti-symmetric solutions, past 
the critical point for the symmetry breaking bifurcation, an asymmetric solution is expected to exist which possesses 
non- vanishing po and pi (as well as a non-zero relative phase between them) , which can be computed from Eq. (j2ip . 
It is anticipated that the presence of loss and gain will not (generically) modify the nature of the bifurcations in 
comparison to the Hamiltonian case (43 |. Namely, an asymmetric solution will bifurcate from the symmetric one 
in the focusing nonlinearity case of s = —1, due to a non- vanishing contribution of the anti-symmetric part in the 
solution, while on the contrary, an asymmetric mode will emanate from the anti-symmetric one in the defocusing 
nonlinearity setting of s = 1 (due to a symmetric contribution within the solution). These results are detailed for a 
particular case example potential in what follows and compared to full numerical results. 



III. NUMERICAL RESULTS 



In our theoretical approximations, the double-well potential is constructed by placing a localized barrier at the 
center of the parabolic trap potential of strength Q,. Particularly, the double-well potential is assumed to be of the 
form: 

V{x) = -n^x^ + Vo sech (-) , (24) 
2 \w/ 

where w is the width of the barrier and Vq its height. The results presented below are for the potential parameters 
fl^ = 0.1, Vq = 5, and w ~ 0.2; we have checked that other parameter values lead to quantitatively similar results . 
For the gain we consider a strength a = 0.2 and a spot size of Xm = 2.0. The damping parameter a is used to vary 
the number of atoms, N, in order to do the continuation. For the above double- well potential, the values of the linear 
eigen-energies are ujq = 0.515729 and uji = 0.677697. The potential setting under consideration is depicted in Fig. [1] 
We have performed a continuation of symmetric, anti-symmetric and asymmetric states in both cases of repulsive 
and attractive interactions. The continuations have been performed by increasing the damping parameter a, which 
is tantamount to decreasing the norm or chemical potential. It is important to note that the chemical potential is no 
longer a free parameter in the present setting in sharp contrast to what is the case in the Hamiltonian regime of atomic 
BECs (see also the discussion of Refs. [13, HO])- Similar results can be obtained by decreasing the pumping parameter 
a. However, a crucial realization that emerges from considering variations of the different parameters is that the spot 
size Xm must be chosen in a very limited range in order for the three above mentioned nonlinear modes to co-exist 
and be potentially stable; outside this range, instabilities lead to breathing multi-bump coherent structures. In what 
follows, the values of Xm = 2 and a = 0.2 have been used unless explicitly indicated otherwise. 



A. Repulsive case 



We start by considering the case of the repulsive interaction with s — 4-1 (and vary a as mentioned above). The 
family of symmetric solutions is found to be always stable. As expected, on the other hand, and in agreement to 
our expectation from the realm of atomic BECs, the anti-symmetric solutions are exponentially unstable for small 
cr, which is tantamount to large polariton population numbers N . They become stable after the symmetry-breaking 
pitchfork bifurcation occurring at a = 1.045 (i.e., for p < pcr = 0.7574 and for N < N^r = 0.5333). The asymmetric 
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FIG. 1: (Color online) The parabolic trapping potential and the localized barrier creating the double- well potential configuration. 
The parameter values used are: = 0.1, w = 0.2, and Vb = 5; the shaded area corresponds to the region where the pumping 
acts, i.e., \x\ < Xm = 2. 




FIG. 2: (Color online) Bifurcation diagrams for the symmetric, anti-symmetric and asymmetric branches for defocusing (re- 
pulsive) nonlinearity (s = 1). Left: Dependence of the chemical potential on the damping parameter. Right: Dependence of 
the (normalized) number of polaritons on the chemical potential. Unstable solutions are depicted by dashed lines on the left 
panel. On the right solid lines display numerical results obtained by a nonlinear (Newton-Raphson) solver of the steady state 
equations of the model of Eq. (IT|, while dashed lines display analytical results of our Galerkin approach. The linear modes are 
located at = 0.5157 and = 0.6667. 



branch that emerges through this bifurcation is stable for /i < 0.7603 and N < 0.5509, i.e., for a narrow parametric 
interval past the bifurcation critical point. However, past this secondary critical point, the asymmetric solutions are 
prone towards an oscillatory instability emerging through a Hopf bifurcation (the critical loss strength in this case 
is (T = 0.989). The relevant bifurcation diagrams are presented in Fig. [5J which shows the dependence of /z on cr, 
as well as the dependence of N on n (note that the latter form of the bifurcation diagram is more commonly used 
in relevant studies). The latter graph also contains the results of the theoretical analysis for the symmetric branch 
of Eq. (|20p and for the anti-symmetric one of Eq. ([T5)) , as well as for the asymmetric branch which is theoretically 
predicted for the parameters of our double- well potential to bifurcate from the anti-symmetric solution for /i > 0.7722 
and N > 0.6661. As can be seen (also from Fig. [5]), there is good agreement between theoretical predictions and 
numerical findings. 

Some case examples of solution profiles for the different branches, together with the results of their corresponding 
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FIG. 3: (Color online) (Left) Real and imaginary part of the wavefunction profile for a symmetric (top) and anti-symmetric 
(bottom) solution. (Right) Their corresponding stability eigenvalues. In all cases a — 1 and the interactions are repulsive 
(s = +l). 
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FIG. 4: (Color online) (Left) Real and imaginary part of the wavefunction profile for an asymmetric solution with a = 0.5. 
(Right) Their corresponding stability eigenvalues. All cases correspond to the repulsive interaction case (s = -|-1). 
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FIG. 5: (Color online) Dependence of the imaginary part of the stability eigenvalues with respect to a for the anti-symmetric 
(left) and asymmetric solutions (right) in the repulsive interaction case (s — +1). 
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FIG. 6: (Color online) Top: Evolution of a perturbed anti-symmetric soliton for a — 1 (left) and a — 0.8 (right) in the repulsive 
case (s = -1-1). The former case relaxes to the asymmetric stationary state, while the latter to the symmetric ground state. 
Bottom: Respective time series for the density at the bottom of the left (solid blue line) and right (dashed red line) wells. 



linear stability analysis as performed by means of the Bogolyubov-de Gennes (BdG) ansatz are shown in Figs. [3] 
andlH The BdG analysis is represented by the speetral plane of the linearization eigenfrequencies uj = Re(cL)) + ilm(a;) . 
Contrary to what is the ease in the Hamiltonian setting of Ref. (where the speetrum is chiefly on the imaginary 
axis), here the spectrum contains predominantly decaying modes with Im(w) < 0. For the stable symmetric ground 
state in Fig. [3l all modes are decaying except for the symmetry mode associated with uj = 0, while for the unstable 
anti-symmetric mode of the bottom panel the eigenfrequency associated with the growth is purely imaginary with 
Im(aj) > 0. On the other hand, for the asymmetric modes of Fig. HI it is evident that shortly past the critical point 
for their emergence, a genuine (now that the system is dissipative, in nature) Hopf bifurcation arises through the 
crossing of a complex conjugate pair through the axis of Im(aj) = 0. Additional Hopf bifurcations happen for smaller 
values of a (larger values of A^), a case example of which is evident in the bottom panel of Fig. SI The dependence 
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FIG. 7: (Color online) Top: Evolution of a perturbed asymmetric soliton with a = 0.8 (left) and a = 0.98 (right) in the 
repulsive case (s = +1). In the former case, perturbation leads to the symmetric ground state attractor, while in the latter 
case, it relaxes to a non-stationary (quasi-periodic) solution. Bottom; Respective time series for the density at the bottom of 
the left (solid blue line) and right (dashed red line) wells. 




FIG. 8: (Color online) Bifurcation diagrams for the symmetric, anti-symmetric and asymmetric branches for focusing (attrac- 
tive) nonlinearity (s = — 1). Left: Dependence of the chemical potential on the damping parameter. Right: Dependence of the 
(normalized) number of polaritons on the chemical potential. Unstable solutions are depicted by dashed lines on the left panel. 
On the right panel solid lines display numerical (Newton-Raphson) results, while dashed lines display analytical (Galerkin) 
results. 



of the imaginary part of the relevant eigenvalues for the anti-symmetric and asymmetric solutions with respect to a 
is shown in Fig. [5J illustrating, respectively, the relevant pitchfork (left panel) and multiple Hopf bifurcations (right 
panel). Naturally, the Hopf bifurcation of the asymmetric branch is anticipated to give rise to a limit cycle attractor 
within the dynamics [the relevant solution is expected to be periodic in the squared modulus of the wavefunction, 
hence quasi-periodic in the original field u{x,t)]. 
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FIG. 9: (Color online) (Left) Real and imaginary part of the wavefunction profile for a symmetric (top) and anti-symmetric 
(bottom) solution. (Right) Their corresponding stability eigenvalues. In all cases a = 1 and the nonlinearity is attractive 
(s = -l). 



Two examples of the dynamics of unstable anti-symmetric solutions are illustrated in Fig. [6l It is observed that the 
unstable solutions generically tend to the stable attractors. However, interestingly, in the a = 1 case, the attractor 
of relevance consists of an asymmetric steady state, while in the cr = 0.8 case it consists of a symmetric one (the 
ground state of the system) . The symmetry and asymmetry of the configurations can be easily seen from the time 
series of the densities and measured, respectively, at the bottom of the left and right wells. These time series 
are depicted in the lower panels of the figure. The relevance of the asymmetric attractor, especially for larger values 
of N (smaller values of a, where the only stable steady state is the symmetric one) is confirmed by the simulation 
shown in the left panel of Fig. [71 where the dynamics of an unstable asymmetric solution is traced, leading indeed 
to the same attractor. The right panel of Fig. [7] shows the evolution of a perturbed asymmetric state close to the 
Hopf bifurcation; in that case, it is observed that the soliton relaxes to a quasi-periodic asymmetric solution. [Recall 
that these solutions have a quasi-periodic evolution for the wavefunction (due to the periodic evolution of the phase 
through e~*^*) but the evolution in density is periodic as the panels show]. That quasi-periodic branch is observed to 
exist in the range a e [0.970, 0.990]. 

B. Attractive case 

In the case of attractive interactions (s = —1), the scenario is similar in nature, except for the origin of the symmetry 
breaking bifurcation. More specifically, now, the asymmetric solutions, which stabilizes at cr = 0.923 = 0.4182 and 
N = 0.7199), bifurcates from the symmetric solutions branch at a = 1.118 {fi — 0.4101 and N = 0.7741). Figures [5] 
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FIG. 10: (Color online) (Left) Real and imaginary part of the wavefunction profile for an asymmetric solution with a = 1 
(top) and a = 0.7 (bottom). (Right) Their corresponding stability eigenvalues. Notice the Hopf bifurcations and the associated 
oscillatory instabilities through two complex pairs which have occurred in the latter case. Here, the nonlinearity is attractive 
is = -l). 




FIG. 11: (Color online) Dependence of the imaginary part of the stability eigenvalues with respect to a for symmetric (left) 
and asymmetric solutions (right). Here, again, the nonlinearity is attractive (s = —1)- 



12 





FIG. 12: (Color online) Top: Evolution of perturbed symmetric solitons with cr = 1 (left) and a = 0.82 (right) for attractive 
nonlinearity (s = —1). Bottom: Respective time series for the density at the bottom of the left (solid blue line) and right 
(dashed red line) wells. 
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FIG. 13: (Color online) Top: Dynamical evolution of the density of the non-stationary asymmetric solution branch found for 
attractive nonlinearity (s = —1) in two cases: a = 0.9 (left) and a — 0.8 (right). Bottom: Respective time series for the density 
at the bottom of the left (solid blue line) and right (dashed red line) wells. 
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and [TT] are the equivalent to Figs. [2] and [5l respectively, but for s = —1. Nevertheless, we observe that both the 
dependence of the chemical potential fi on the nonlinear saturation parameter a and that of N on /x is, in fact, 
non-monotonic for this example in the case of the bifurcating asymmetric branch. This clearly indicates (sec the 
right panel of Fig. [5]) that the relevant bifurcation is subcritical (as the chemical potential /i is decreased, which is the 
natural direction of variation off of the linear limit). This is contrary to the corresponding supercritical expectation of 
its Hamiltonian analog [i^jIElj. It should be noticed, however, that other examples where such subcritical bifurcations 
have been previously reported in Refs. [ssl . [s^ although in neither case was the nonlinearity purely cubic as was the 
case here (and they did not contain driving/damping effects). Importantly, it should also be pointed out that the 
analytical prediction of the Galcrkin approach suggests a supercritical scenario for fi < 0.4247 and N > 0.6590. 
Despite the inability of the approximation to capture the short subcritical segment of the bifurcating branch, we 
nevertheless see that the Galerkin method is a useful tool for obtaining an estimate of the relevant critical point. 

An additional feature worth pointing out concerns the nature of the instabilities of the different branches as detailed 
in Figs. |9] and [TOl While the symmetric branch becomes unstable at the relevant critical point by developing an 
imaginary eigenfrequency with Im(w) > (the rest of the spectrum has Im(cj) < 0), the anti-symmetric state remains 
dynamically robust. On the other hand, the asymmetric branch emerges as stable at the critical point of the symmetry 
breaking but shortly thereafter (for a < 0.923), it becomes subject to a Hopf bifurcation through the crossing of the 
axis with Im(ci;) = of a complex eigenvalue pair. In fact, for a < 0.74, a secondary Hopf bifurcation has occurred 
and is mirrored in the two complex pairs with Im(a;) > shown in Fig. 1101 This phenomenology is enforced by 
Fig. [TT] which illustrates the dependence of the relevant stability eigenvalues on the nonlinear loss parameter a (see 
the right panel for the sequence of Hopf bifurcations, while the left panel highlights the symmetry-breaking induced 
crossing of a single eigenfrequency pair for the symmetric branch). As in the repulsive case, the Hopf bifurcation of 
the asymmetric branch is anticipated to give rise to a limit cycle attractor within the dynamics. 

The dynamics of Figs. [1^ and [T2] naturally reflects the above conclusions. In particular, the evolution of the 
symmetric state in the double-well potential of the left panel of Fig. [12] gives rise to the asymmetric state as the 
latter is stable and indeed an attractor for the value of cr = 1. The right panel of the figure displays the evolution of 
a perturbed symmetric solution tending to an anti-symmetric one; in that case, the asymmetric solution is unstable 
and no longer a dynamical attractor. 

On the other hand, Fig. [13] shows different case examples of the (unstable via the Hopf) asymmetric branch for 
different values of cr. In those cases, the asymmetric branch is no longer a stable stationary state and as a result the 
dynamics becomes periodic in the modulus (quasi-periodic in the original field) for a S [0.74,0.92]. It is interesting 
to follow the changes in the dynamics for these periodic states as a is decreased below the bifurcating point from 
the asymmetric branch. In particular, close to bifurcation point, the periodic evolution remains proximal to the state 
from which it emanates, namely the asymmetric state as it can be seen in the left panels of Fig. 13. However, as a 
is decreased further from the bifurcation point, the instability of the asymmetric state is stronger and the departure 
from the asymmetric solution is more significant. In particular, it is interesting to notice that for smaller values of cr, 
the solution tends to display strong oscillations of the densities resembling the symmetric tunneling of matter from 
one well to the other. An example of this evolution for a = 0.8 is depicted in the right panels of Fig. 13 where it 
is evident that the oscillations in the two wells become similar to each other but with a phase shift between them, 
leading to an effective rc-symmetrization of the dynamics. 

It is also interesting to highlight here the difference between the repulsive case of Figs. [6] and [7] and the attractive 
case of Figs. [T^ and \W\ In the former case, when the emerging asymmetric branch is unstable the dynamics typically 
is found to lead to the stable ground state of the system (the symmetric one). On the other hand, for the attractive 
case, when both the symmetric and the asymmetric branch are destabilized, the dynamics does not resort to the 
excited (yet stable) anti-symmetric state. Instead, it leads to periodic oscillations in the density between the two 
wells. 

Finally, we have considered the effect of varying the spot size fixing a — I. In the repulsive case, the symmetric 
branch is stable for Xm G [0.9,5.7]; out of this range, the instabilities are caused by a Hopf bifurcation cascade 
and develop into non-stationary multi-dark soliton waveforms, similar to the states that were previously reported 
in Ref. (soj (but for a purely parabolic trap). The anti-symmetric branch, which is unstable for every Xm (for this 
value of cr), experiences a bifurcation cascade for Xm < 2.0 and Xm > 5.3. The instabilities at Xm G (2.0,5.3) are 
the exponential ones previously explored. However, considering higher values of cr, a stability range appears which is 
enlarged for growing cr. A similar effect is observed for the asymmetric branch, i.e., there is a small stability interval 
Xm S [1.9,2.0] that is enlarged when cr is decreased. Outside this range, the branch experiences Hopf bifurcation 
cascades. 

The above mentioned scenario is almost equivalent for the attractive case, except that the symmetric and anti- 
symmetric branches are interchanged. In that case, the anti-symmetric branch is stable for Xm G [2.0,4.8]; the 
symmetric branch is now stable for Xm G [1-0, 1.9], starting the Hopf cascade at Xm = 4.5. The asymmetric branch is 
stable for x„i G [1.0, 2.0], while being oscillatorily unstable for other values of Xm- 
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IV. CONCLUSIONS AND FUTURE CHALLENGES 

111 the present work, we studied the existenee of solutions, their speetral stabihty and nonhnear dynamics for the 
case of a polariton condensate confined in a quasi-lD double well potential. Motivated by recent developments for 
the study of polaritons in such settings (23l - |29l . |47| , and by the work of Ref . [3 which proposed a two- well model, we 
presented a systematic Galerkin analysis for the model with the gain over a localized spot and nonlinear saturation 
loss formulated in Refs. It was theoretically predicted that nonlinear states emanate from the corresponding 

linear ones of the potential and that bifurcations are expected to arise, similarly to the Hamiltonian analog of this 
setting studied earlier in the context of atomic BECs. Such symmetry breaking pitchfork events emerge from the 
anti-symmetric, first excited state in the case of the repulsive interactions, while they arise from the symmetric ground 
state branch in the case of attractive ones. Despite the similarities with the atomic BEG case, nontrivial differences 
exist as well. One of them concerns the nature of the bifurcation, which in the attractive case was found to be 
weakly subcritical (instead of supercritical) upon decrease of the chemical potential. Importantly also, the resulting 
asymmetric branches aside from narrow intervals of stability are generically found to be unstable due to genuine Hopf 
bifurcations, which, in turn, give rise to periodic orbits (in the density). While in the repulsive case, the dynamics of 
anti-symmetric and asymmetric branches is found to be attracted to the ground state when both of them are unstable, 
the periodic orbits are essential to the evolution in the case of attractive interactions as they seem to constitute the 
robust dynamical attractor. 

This is merely the first step in the examination of the similarities (but also the differences) of the polariton BECs 
and their atomic counterparts within a setting that contains the interplay of a double-well potential and nonlinear 
interactions. Yet, our study paves the way for a number of potential future avenues. On the one hand, one can consider 
the more detailed model of Refs. p34l6| and examine whether the inclusion of the diffusive dynamics of the exciton 
population induces any qualitative differences in the features reported herein. On the other hand, and bearing in mind 
the predominantly two-dimensional nature of the polariton dynamics, one can envision generalizations of the potential 
considered herein in a 2D realm. Relevant possibilities may include not only the straightforward generalization of 
a double well encompassing two quasi-one-dimensional tracks, but also that of a genuinely two-dimensional four 
well potential that has recently been examined in detail in atomic BECs [s^l- Even in the context of the present 
model, there are further possibilities to explore, including the systematic investigation of the emergent periodic orbits 
and their Floquet spectral stability analysis. Such studies are currently in progress and will be reported in future 
publications. 
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